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GUST RESPONSE OF COMMERCIAL JET AIRCRAFT 


INCLUDING EFFECTS OF AUTOPILOT OPERATION 

By Joseph H. Goldberg 
Villanova University 

SUMMARY 

This paper presents a simplified theory of the vertical c.g. 
acceleration gust response of aircraft. The development is based 
on a three-degree of freedom model including pitch, vertical dis- 
placement, and longitudinal control deflection resulting from the 
autopilot pitch and vertical displacement feedbacks. Higher-order 
autopilot transfer functions are utilized to provide the necessary 
accuracy in the prediction of gust response characteristics such 
as the gust response parameters relating gust input to accelera- 
tion, the acceleration zero-level crossing number, and power spec- 
tral density properties. The method is applied to four representa- 
tive commercial jet aircraft for a broad range of operating condi- 
tions and a comparison of all gust response characteristics is 
made. It is shown that autopilot operation relative to the controls 
- fixed case causes response attenuation of from 10 percent to 
approximately 25 percent depending on flight condition and increases 
in crossing number up to 30 percent, with variations between air- 
craft of from 5 percent to 10 percent, in general, which reflect 
the differences in autopilot design. A computer program descrip- 
tion and listing of the calculation procedure suitable for the 
general application of the theory to any airplane - autopilot 





INTRODUCTION 


combination is presented. 

Considerable effort, by industry and government agencies, 
has been devoted to the gust response problem Iqieginning with the 
first NACA report. Currently the Safety and Operating Problems 
Branch at NASA Langley is engaged in an extensive study of digital 
flight recorder data from commercial jet operations. Central to 
this effort is the analysis of vertical c.g. acceleration records 
for the purpose of isolating the gust velocity input. A two- 
degree-of-freedom theory, unique in approach (ref. 1) , was employed 
for the reduction and analysis of these data. (Additional impor- 
tant work is contained in references 2 through 9.) It was found 
(ref. 10) that a three-degree-of- freedom theory was required to 
account accurately for the influence of autopilot operation on the 
airplane gust response. 

The study was undertaken for the purpose of extending the 
theory of reference 1 to include the longitudinal control auto- 
pilot; and to develop working methods for the prediction of gust 
response characteristics such as the gust response parameter A^, 
zero-crossing number Nq, and the power spectral density PSD, for 
various modes of autopilot operation. Four representative commer- 
cial jet designs were investigated over a wide range of flight 
conditions, loadings, and autopilot operating modes. The auto- 
pilots were modeled in a high-order linearized version to produce 
complex interactions with the airframe dynamics if present. The 
theoretical development included elastic airframe aerodynamics 
including compressibility effects over the 
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frequency range of interest. A computer program (appendix) was 
developed specifically for this purpose and structured for effi- 
cient calculation of the gust response parameters for each design, 
with varying aerodynamic and autopilot inputs corresponding to 
different combinations of gross weight, Mach number, altitude, and 
autopilot mode. 

SYMBOLS 

a wing lift-curve slope, rad“^, based on wing area and ‘wing 
local flow direction; real nxamber coefficient in the air- 
plane y transfer function 

a imaginary number coefficient in the airplane y transfer 

function 

Sq complete airplane lift-curve slope, rad“^, based on wing 
a,rea and wing local flow direction 
at horizontal tail lift-curve slope, rad“l, based on horizontal 
tail area and tail local flow direction 
A real number coefficient in the autopilot transfer function; 

subscripts indicate independent variable 
Ag ratio of normal acceleration increment to continuous gust 
velocity, g units/m/s (g units/fps) . 

A^ wing geometric aspect ratio 

At tail geometric aspect ratio 

A/P autopilot 

APF autopilot factor, for A(j or Nq 

b real number coefficient in the airplane y transfer function 
b imaginary niimber coefficient in the airplane y transfer 

function 
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B imaginary number coefficient in the autopilot transfer 

function; subscripts indicate independent variable 
c real number coefficient in the airplane y transfer function 
c imaginary number coefficient in the airplane y transfer 

function 

c wing mean geometric chord, m (ft) 

horizontal tail mean geometric chord, m (ft) 
e effective wing moment arm, or distance from wing lift 

aerodynamic center to c.g. as required for S.M. simulation, 
positive forward of c.g., m (ft) 

e-t horizontal tail distance from tail aerodynamic center to c.g. 

positive aft of c.g., m (ft) 
f frequency, Hz 

f^ generalized airplane transfer function, 

fc upper integration limit, Hz 

= 3 at/8 6 

tail 3at/95 factor 

g acceleration due to gravity, 9.81m/sec2 (32.2 ft/sec^) 

hp pressure altitude, m (ft) 

i imaginary number, 

k nondimensional gust frequency, CO c/2V 

kc nondimensional upper integration limit, TT cfc/V 

kg nondimensional acceleration zero-level crossing mamber, 

Tf cNq/V 

j(_ pitch radius of gyration, m (ft) 

K modified gust alleviation factor, based on rather than 

®w 
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k' K integral weighted by KkQ 

Kj 4 compressibility correction factors for Aa, single or double 
primed 

L Von Karman tubulence scale m (ft) 

Li total wing lift, at c/4, N (Ibf) 

L 2 total horizontal tail lift, atc^./4, N (Ibf) 
m airplane total mass, kg (slugs) 

M Mach number 

n normal, i.e. vertical acceleration, g units 

An normal acceleration increment, n-1, g units 

Ng niimber of positive slope zero-level crossings of normal 
acceleration increment per second 
PSD power spectral density 

rg wing lift slope ratio, a/ao 

weighted horizontal tail lift slope ratio, St^t/Sag 

r 2 horizontal tail nondimensional downwash point distance, 
2et/c + c^/c + i^(2e/c-l) 
rms root mean square 

s Laplace complex variable, ioj 

Sj^ nondimensional wi distance from W 2 , 2e/c~l + 2et/c + ct/c 
S wing area, m^ (ft^) 

S.M. airplane controls - fixed static margin, fraction of c 
horizontal tail area, m2 (ft2) 

u input gust velocity amplitude, positive up, m/s (ft/sec) 

Ua continuous gust velocity, m/s (ft/sec) 

V true airspeed, m/s (ft/sec) 
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wing downwash, at 3c/4, positive up, m/s (ft/sec) 
horizontal tail downwash, at 30^/4, positive up, m/s (ft/sec) 
vertical displacement (inertial frame) of the c.g.,m (ft) 
vertical displacement of wing downwash point, m (ft) 
vertical displacement of tail downwash point, m (ft) 
wing angle of attack, measured from local flow direction 
to zero lift line, rad 

tail angle of attack, measured from local flow direction 

to zero lift line, rad 

Glauret compressibility factor , l-M^ 

longitudinal control deflection, positive trailing edge 
down , rad 

gust penetration time, sec 

tail downwash derivative, 9 e/a a, where e is the downwash 
angle 

mass ratio parameter, 2m/aQpSo 
ambient air density, kg/m^ (slugs/ft^) 
rms acceleration, g units 
rms gust velocity, m/s (ft/sec) 

rms modified gust velocity. O’ (2L/c:) ^ , m/s (ft/sec) 

airplane pitch angle, corresponds to wing angle of attack, 
rad 

gust PSD, m?42/Hz (ft^sec2/Hz) 
acceleration PSD, (g units) 2 /hz 
gust input frequency, rad/sec 



EQUATIONS OF MOTION 


The equations of motion are for a two-degree of freedom 
body with longitudinal control (elevator or full stabilizer 
motion as the primary contrdl) . A model based on two simplified 
lifting surfaces for the wing and horizontal tail is adopted. 

The analysis includes first order nonsteady lift effects, over- 
all corrections for compressibility, first order corrections 
for gust penetration and gust lag effects on the tail relative to 
the wing. Since the primary frequency range of interest is low, 
downwash effects are included on the basis- of a quasi-steady state 
condition.. Static aeroelastic effects, are considered by using 
aerodynajnic characteristics at the equilibrium condition for the 
elastic airframe. 

The equations of motion, provided by the author of reference 
1, are extended iverein to include an autopilot actuated control 
surface with autopilot modeling in a highrorder version. The equa 
tion development follows that of reference 1 .but is simplified by 
the use of a quasi-steady nonstea:dy lift assumption, Downwash 
expressions are given in terms of the threer-degrees of freedom 
obtained when including control deflections. Solutions for ver-^ 
tical acceleration at the. c.g, are developed in a closedform for 
the sinusoidal vertical gust input of unit magnitude.. From these 
frequency response functions, spectral results are derived. 

Referring to figure 1, the vertical displacements at wing 
and tail downwash control points are seen to be given by 

Yl = y + (e - |)(|) 

Yo = y - (e . + <j> 
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The lift forces are 

Li = I /)SVwi 

^2 ~ J)S^VW2 

where the wing downwash is 

wj_ = V<}) - Y]_ + u 
and the tail downwash is 

• • 

W2 = V((|) + F^6) ~ Y 2 + ^5 Y~ ^ ~ T^^l 

In this last expression the factor is equivalent to the elevator 
or control lift effectiveness, i.e., ^“t /96 (unity for an all-move- 
able surface) , and Fg is taken to be zero for elevator control and 
unity for an all-moveable primar;y control surface. Lift curve 
slopes for the wing and tail, a and a^, are for the isolated com- 
ponents based on their individual areas and are the actual three- 
dimensional quantities at the equilibrium operating condition as 
required for this approach. The vertical gust input is u at the 
wing downwash point and is delayed by e time when reaching the tail. 
A major effect, the steady state wing downwash velocity at the tail 
control point is given in terms of denoting the familiar 3 e/a a 
where e is the downwash angle. 

Taking the time derivatives of yj^ and y2 and making substitu- 
tions leads to 

* dr * 

Wj^ = V<|) - y - (e - 2-) <|) + u 

• • • 

W2 = V(<f + F56) - y + (e^ - + (u)t_e 

• 2 2 

-T)[v(t> - y-(e-|-)<f>+v^ 
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Now the vertical force and pitching moment equations are 
developed. The vertical force is given by 
my = Lj + L 2 = f jjSVwi + ^ fS^Vw2 

The pitching moment equation for the complete configuration , 
(including fuselage and mutual interference effects) can be 
expressed in terms of wing and tail by making the wing lift 
moment arm, e, that required to duplicate the moment character- 
istics of the overall configuration. The moment equation is then 


m^^^ = eLi - etL2 

- e| pSVwi - e ^-^t j)StVw 2 


Nondimensionalizing the force and moment equation results in 


2m 


V = ^ '^l + 


W' 


aQ ^0 V ^ S V 

where ap is the complete airplane lift curve slope. 
^2 


— -2m^ — X = ® a ''I _ et 
anPSV^ c ^ c ^ — 


2m 

Also, aQpSV* 


‘0 V 


ao S 


w^ 

V 


= 


where p, the mass ratio parameter^ is 
2m 

V = agpSc 


Thus, 

aQpSV^c v2 


The final simplified forms for the force, moment, wing down- 
wash, and tail downwash in terms of the five variables considered 
at this stage are now derived. Here the equations are written in 


9 



the frequency rather than time domain. Thus the input variable u 
represents the amplitude of a sinusoidal vertical gust of frequency 
CaJ. Also to be noted is that the output or frequency response of 
interest is the vertical displacement y. (The lowest order in y 
in the . preceeding is y, and the magnitude of the vertical acceler- 
ation frequency response is Cj|y|; but the control law for 6 will 
be defined in terms of y and (f>.) The four equations are 


(fi £ iJ)Z + (ro)!li + (r,)!^ = 0 
V V ^ V ■‘■V 

(p ^0)2)4, + - (®tri)'^2 = 0 

v2 . ° V c- ^ 

^y + r (e-c)iW - ll ^ = u 

V ‘-2 V~ A7 Tf 


V V 

c 




(1) 

( 2 ) 

(3) 




- (Fg + j-CJ) 6 + W2 = (e~^^^Tt)li 

V V 


where r 


0 . 


a 

ao 


and r, = ^ 

ao S 


(4) 


Now the forcing frequency is made dimensionless by normalizing on 
2V/C, thus the nondimensional input frequency is k = « 

Also the penetration time interval e is given by 

. = 1 


= ^(e-|+et+^) 

Then e“^^®l = cosks^^-isinks^ 


where si = (^) e = 2®. - 1 + 
c c 


2&J. c 

+ t 


The four equations now reduce to 
|^ 2 pk 2 {—)] y + rQWi + riW 2 = 0 

k2j V(|) + ® rQwi - ^ riW2 = 


( 1 ) 

( 2 ) 
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(3) 


[k(^)i] y + [(^-l)ki-l3 V(j) + wi = u 

l^k(l-r^) y - [](l-n)+r 2 ki] V(|) - + Fg^ k:^V6 

+ W2 = [^(cosksj^-i^) - (sinks^^) ij u (.4) 

where T 2 = + 2Jt + (^-1) 

c c c 

The control law (the equation for V6 in terms of the transfer 

functions for V6 due to y and due to V(J)) is now developed for 

introduction into equation (4) . It is initially given here in the 

Laplace variable s, later to be replaced by ico, or in nondimensional 

form as (^)ki. 
c 

V6 = {Yl)y + (Y±)V^ 
y V(J) 

where the bracketed terms are the two control transfer functions . 

As used in this study they are the linearized (limiter, dead-zone 
and hysteresis operations omitted) functions for the individual 
autopilot designs for each airplane. The general form found 
applicable to higher-order modeling for all autopilots, as given 
in terms of the computer program notation (see subroutine CALC in 
the appendix for the definitions) , is 

V6 = BD [(YN) + 

DN L yd PD J 

Introducing the control effectiveness function, we now define the 
following 

Ay + Byi = -[fs + (Ffift k)i1 (5E) (YN) 

c DN YD 

A. + B.i = -Ff 5 + (rjEt k)i"| (DD) (PN) 

^ ^ L c -I DN PD 
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Equation ( 4 ) , with V6 eliminated, then becomes 

[k(l-T^) (|Y)i + (Ay+Byi)]y 

- U(l-r\) + r2ki - (A<j,+B<j)i)J V4, + W2 

= Q(cosksl-i|) - (sinksi)i^u ( 4 ) 

The preceding set of four equations is readily reduced to a 
2 x 2 matrix by solving ( 3 ) for wi and ( 4 ) for W2 and substituting 
in ( 1 ) and ( 2 ) . With u set to unity for a unit gust input ampli- 
tude, the y solution is found to be 

y = (cj+ci) (b2+b2)- (C2+C2) (bi+bj) 

(ai+ai) (b2+b2)- (a2+a2) (bi+bi) 

where 

ai + ai = {riAy-2y (^)k2} + ( riBy+ {^) j]’ rQ+ri (1-t^)] k)i 

c c 


^2 

+ 

a2 = 

{-et I 

c 

•lAy} + 


f ^0- ^ rid-il)] k}i 

bl 

+ 

bi = 

{riA^ 

- [ro+ri (1-1^)] } + (riB^+ [ro (|^-1) -rir2] k}i 

^2 

+ 


{-et. I 

c 

[ 

e ro - ct £r£l-T\) 
- c c 

+ 4 v ^2 j^2l } 

c2 J 




+ 

c 

L riB<j, 

+ [e ro(^-l) + 
L c c 

®t rir2”lk}i 
^ -J 

Cl 

+ 

Cl = 

[ro + 

ri (cosksi-r|)J - [ri(£ 

iinksi^ i 

C2 

+ 

C2 = 

p ro 
^ c 

- 

c 

ri (cosksi-q^)^ + [ 

^t ri (sinksi)J i 


The y solution obtained here is the frequency response com- 
plex amplitude for a given input nondimensional k frequency. From 
this solution all of the gust response characteristics of interest 
in the study are obtained as described in the following section. 
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CALCULATION OF GUST RESPONSE CHARACTERISTICS 
The vertical velocity frequency response is developed from 
the vertical displacement with 
y = i y = i(j (^)ky 

Expressed in real and imaginary parts, 

• • » 
y = Yr + iyi 


It should be noted that this 

quantity is nondimens ional since y is per unit u, each in distance 

per time units. The absolute vertical velocity magnitude square 

(as required for the generalized transfer function below) is then 
• • • ' • • 

|yp = (Yr + iyi) (yr - iyi) 

An intermediate nondimensional response quantity termed the 
generalized airplane transfer function (developed in reference 1) 
and denoted as f^, is first determined from 

f]_ = 4y 2^2 I y [ 2 

(1 + k + ) (1 + 0.55A^) 

Where = i-m^ and is the wing aspect ratio. In this equation 
the denominator factors provide the necessary corrections for gust 
penetration and spanwise gust distribution effects. 

The generalized gust response PSD is obtained from the product 
of the airplane transfer function and the normalized gust input 
PSD. This gust input PSD is discussed in reference 1, and is 
given as 

4 >w = (^) [1 + |(1.339 ^ k)2] 

^ [1 + (1.339 ^ k)2]ll/16 
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where the equivalent rms gust 0 ]^ is related to the actual rms 
'^w by 

2 _ 

'^l ,2L '2/3 

7T (—) 

The turbulence scale L, a fundamental parameter but one having 
less influence at values of 2L/c in excess of about 100 (ref. 1) , 
is taken as 305 m (1000 ft) throughout this study. 

Response spectra are presented as the product of frequency 
and the PSD (affording an indication of the frequency intervals 
containing the greatest energy content) . This generalized pro- 
duct of frequency and PSD for the acceleration response to gust 
is designated GFPSD and given as 
GFPSD = k(f, ) 

Also, this quantity is identical to that contained in reference 1. 

The equivalent gust alleviation factor is now found. It is 
so termed because it is based on 0 ^ rather than a^. Designated 
by K it is calculated from the following integral. 




K 


where k^,, the upper limit, is the nondime nsional cut-off frequency 
given as 


^fc 

Here fc depends upon the specific airplane, and the value corres- 
ponds to that used in the reduction of flight response data in 
the Langley Digital VGH program. In general it is about five 
times the airplane's controls 
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- fixed short period natural frequency. This criterion permits 
inclusion of all spectral content contributing most significantly 
to the c.g. acceleration response. It should be noted that these 
alleviation factors are considerably greater than those based on 
the actual rms gust a^. The alleviation factor is increased by 
the term 

[. (|£) 2 / 3 ] '/2 

which is approximately 8 for the average configuration. 

The gust response parameter A(j relating the incremental load 
factor to the continuous gust velocity is defined as 



and is determined from 

^ (^)K 
g]i 

The two normalized frequency times PSD quantities of major 
interest are the frequency PSD products based on rms gust velocity 
and based on rms acceleration. The first is called the gust 
normalized frequency times PSD, or GNPPSD, and is obtained from 

GNFPSD = GFPSD 

Vgy / 'ir^2L^^j^Z3 

c 

and can be denoted as 

GNFPSD = • (? 2 ( ( 

The second is called the acceleration normalized frequency times 
PSD, or ANFPSD, and is obtained from 
ANFPSD = 
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and can be denoted as 


ANFPSD = , nondimensional 


The zero level crossing number Nq is the number of positive 
slope crossings per second of the zero-acceleration increment 
time history. It is found from 


No 


= V/c 


ko 


where 


and 


ko = 
u ^ 


K 






In summary, the theoretical development undertaken so far 
has shown the dependence of the response quantities of importance, 
such as Ag, Nq/ GNFPSD and ANFPSD, on the fundamental equation set 
comprising wing downwash, tail downwash, total lift, and total 
pitching moment, in conjunction with the autopilot control law. 
Inputs for the specific configurations to be considered are ex'- 
pressed in generalized and nondimensional form. Autopilot inputs 
are defined by gains. and time constants in linearized transfer 
functions. These are kept in a higher-order version for the 
modeling process to retain the effects of all possible airplane- 
autopilot interactions. 


AIRPLANE CHARACTERISTICS 

Airplane inputs can be classified in three distinct groups; 
those constant for all cases for a given design; those that apply 
to one design only but that vary with flight condition (altitude 
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and Mach number) ; and those that are a function of the nondimen- 
sional forcing frequency k as well as flight condition for a 
particular airplane. Common to all configurations and considered 
constant throughout is the c.g. location at 0.25c and the turbu- 
lence scale of 305 m (1000 ft) . Also, a standard atmosphere has 
been assumed. The airplane cut-off frequency f^ varied from 0.75 
Hz to 1.50 Hz depending upon the configuration. 

The complete airplane lift curve slope ag is used in combi- 
nation with the wing aspect ratio, tail area ratio, tail aspect 
ratio and downwash derivative to calculate the isolated component 
wing and tail lift curve slopes. The two equations so employed 
are 

ag = a + at ^(.l-'i|) 
at " a 


1 + " ^ {-i -^ji) 

0 , 9it "t A^ 

In this manner, consistency is maintained between all three slopes, 
which, of course, figure independently in the analysis. Also, the 
wing and tail contributions to the complete configuration static 
longitudinal stability is used to determine the effective wing 
lift arm e from the following equation 


- =1^1^ <^-a) (s.M.) 


c ° fc) c 

where S.M. is the controls-f ixed static margin, in fraction of 
chord, for the specific airplane and flight condition. The 
controls-f ixed short period mode natural frequency is found, in 
nondimensional form, from 



ksp = T2^ = 27tW [l-2(^)(|£H|t)2+u(S.M.)]'i 

In the preceding, hq, and S.M. are considered to vary with 
flight condition and are taken at the corresponding altitude and 
Mach number for the elastic airframe. 

Flight conditions for a given airplane are defined by pressure 
altitude, Mach number and gross weight. For each design a number 
of altitude and M combinations (varying from 7 to 14 depending on 
the individual case) were selected based on complete (take-off to 
maximum cruise to landing) flight profiles for each. Additionally, 
3' gross weights were chosen for each condition as that correspond- 
ing to minimum, average, and maximum operational values. This 
approach lead to an excessively broad range of p values for the 
analysis when considering the more restrictive conditions for 
autopilot operation but it was thought to be advantageous for 
correlation purposes. 

A range of altitudes and Mach numbers was specified in order 
to identify three distinct operating conditions which would be 
applicable to all four airplanes. These are termed and defined as: 
Climb/Descent, M from 0,50 to 0.58, altitude from 3,050 m (10,000 
ft) to 6,100 m (20,000 ft); Low Cruise, M from 0.80 to 0.81, 
altitude from 8,530 m (28,000 ft) to 9,140 m (30,000 ft); High 
Cruise, M from 0.84 to 0.89, altitude from 9,750 m (32,000 ft) to 
12,220 m (40,000 ft). Table 1 presents the nondimensional geo- 
metric and aerodynamic parameters required for each airplane in 
the study. Noted are the nonvarying quantities ^t/g, "^t/c, ®t/c 
and ^/c . At each general condition, the nondimensional wing and 
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tail slopes, ^/ag and ^t/ag, are given as well as the tail lift 
effectiveness ri. Also included is the nondimensional effective 
wing moment arm ®/c as required to simulate the overall static 
stability. Lastly, the nondimensional cut-off frequency kc, as. 
defined in the integral limits, is given. 

The autopilot modes of operation considered for the study 
are those which are appropriate to each design for climb, cruise, 
turbulent setting in cruise, and descent. It was found that all 
of these incorporated pitch-hold (with altitude-hold in cruise) 
with gains and signal shaping or filtering as required. Climb 
and descent involved altitude-rate as well in some cases, and two 
had provision for reduced gains in the turbulent setting. Three 
standard autopilot modes are then adopted and designated as: a 

combined climb and descent condition (Climb/Descent) ; basic cruise 
with pitch and altitude-hold (Cruise); and cruise with reduced 
gains excluding altitude hold (Turbulent) , The autopilot-alone 
characteristics for each of the airplanes is illustrated in fre- 
quency plots for each autopilot. These have been evaluated as the 
log magnitude of the effective change in tail angle of attack 
(due to control deflection) with pitch angle and with vertical 
displacement, and are given by 

^ l°5ig IA 4 , + B,j,i| 

and 

I (.k). , , 

l-y^l - 20 log^Q I Ay + Byil 

where the A and B coefficients are as previously defined. It 
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should be noted that these quantities are not referenced to zero- 
frequency levels, hence they are not amplitude ratios in the usual 
sense. To permit a more general comparison of the individual auto- 
pilots, including the importance of static loop gains, they are 
so defined, as well as by including the aerodynamic control effect- 
iveness. The frequency response results are presented in figure 2 
as. a function of nondiraensional frequency (ratio of forcing fre- 
quency to the short period frequency) for the Low Cruise flight 
condition in the Cruise autopilot mode. 

The frequency responses shown for each autopilot transfer 
function (pitch feedback and altitude feedback) provide a means of 
evaluating basic autopilot differences. Immediately evident is 
the fairly close correspondence between all pitch channels in the 
frequency band spanning the short period natural frequency. Larger 
differences are noted between the altitude channel frequency res- 
ponses, but in general these are only secondary to the pitch fre- 
quency responses. It should be appreciated that only the combina- 
tion of airframe and autopilot can determine the particular air- 
plane gust characterization. 

GUST RESPONSE PARAMETER Act 

Initial studies were made of the A^ parameter as a function 
Of the principal operational variables. It was determined that 
regression functions were obtainable with about 1 percent average 
error when the independent variables were V/c, m and y. Data 
correlations for each airplane were performed using equations of 
the form 
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gAg ^ Cl 

km(v/c) 

where the compressibility correction Kf^ is closely fitted with 
Km = _£2_ 

32 

and Cl, C2r and C 3 and n are constant for a particular design and 
autopilot mode combination. Significant here is the fact that 
these apply over the broad spectrum of flight encompassing take- 
off, climb, cruise, descent and landing at all gross weights; and 
also that the method is applicable to each of four autopilot modes . 
The procedure used was to generate first these regression functions 
for each autopilot mode and then derive an average compressibility 
correction factor for all airplanes at a given mode. The results 
so determined are presented as a function of y in figures 3 through 
6 , in which the average compressibility correction is denoted as 
Km^ . Examination of the results here exhibit a close correspondence 
of the Ag parameter for the different designs, with the exception 
of the low y region in the Cruise mode. Since this is the one 
autopilot mode whose design would likely be optomized for cruise 
altitudes of 7,620 m (25,000 ft) or higher, autopilot effects at 
y near landing may differ significantly. 

Figure 3 for the fixed-control case can't be compared directly 
to the data in figures 4 through 6 for a direct indication of Aa 
attenuation by the autopilots since each has a different Km^ 
function. However, the attenuation effects are quantitatively 
perceived when an autopilot factor APF is defined, as the ratio of 
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A (j with A/P on to that with A/P off. The APF, denoted A^ APF, is 
determined by the ratio of the ordinate value in either figure 4 , 
5 or 6 to the ordinate in figure 3 at the same y . This approach 
results in the determination of the parameter 




APF 


f (y) 


and the ratio of compressibility factors, denoted is in 

general of the form 



C4 - 

Cg - m2 


where C 4 , C 5 and Cg are additional constants for a given autopilot 
mode. The APF normalized now by the Kjyj^ compressibility factors 
are presented in figures 7, 8 and 9. Here the APF is shown in two 
different representations - the APF for any M and the APF for the 
specific M appropriate to that autopilot mode. The curves indi- 
cate that attenuation of gust response, increasing with y, is 
achieved throughout except for one airplane at the extreme low end 
of the y range in the A/P Cruise mode. Additionally, it is evi- 
dent that a general consistency of trends exists for the APF para- 
meter of the different designs considering the distinctive nature 
of the autopilot frequency responses as previously shown. 

In summary, over the operative y range for each aircraft, Ag 
reductions due to autopilot and as a function of y are; between 
4 percent and 30 percent in the Climb/Descent mode; between -13 
percent and 36 percent in the Cruise mode; and between 9 percent 
and 31 percent in the Turbulent mode (based on two of the four 
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aircraft) . 


ZERO-LEVEL CROSSING NUMBER Nq 

An initial investigation of all Nq results was conducted in 
order to establish the basic dependencies of Nq on the operational 
variables, as was done for Aa . It was found from regression ana- 
lysis that V/c and y (M omitted) were the governing quantities but 
in a more complex manner than for the Ao parameter. These fitting 
functions, semi-logarithmic in V/c and y, were employed to generate 
curves of Nq versus y at selected ^/o values for each of the A/P 
modes, Nq results at V/c of 20 sec”l, 30 sec“l, and 40 sec“^ are 
presented in figures 10, 11 and 12 with the A/P off. These V/c 
magnitudes span the operating range from climb to high Mach number 
cruise conditions. It is seen that Nq increases with V/c and de- 
creases with y, with an overall Nq variation of from 0.2 to 0.6 
crossings per sec. The variation of the Nq level among airplane 
designs for e given condition is greater than it was in the case 
of Aa , as might be expected on the basis of the two-parameter (V/c 
and y) versus three^parameter (V/c, M and y) correlation.. 

With regard to the autopilot effect on Nq, again an APF is 
employed representing the ratio of Nq with A/P on to that for A/P 
off (denoted Nq APF) at the same conditions of V/c , y and A/P 
mode. These factors however were not obtained from generalized 
Nq values with A/P on (as with Aa) but rather were derived by a 
regression analysis of APF directly from the computer output. 

This procedure proved to be more convenient considering the added 
complexity of the Nq fitting functions. Results so derived are 
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shown in figures 13 through 18. The Climb/Descent A/P mode 
factors are given in figures 13 and 14 for ^/c of. 20 sec“l. and 
30 sec~l, respectively. Except for one airplane in a limited low 
p region, Nq APF exceeds unity, as expected, with Nq increases 
indicated of 10 percent to 32 percent depending on airplane con- 
figuration, p, and ^/c. An increase in APF with p is evident, and 
the V/c effect is not well defined but is generally limited to 
only a few. percent at a given p. In the Cruise A/P mode results 
are presented for ^/o levels of 30 sec“l and 40 sec“l in figures 
15 and 16, respectively. In this autopilot mode wider variations 
in Nq APF are shown with Ng increases of from zero to 37 percent, 
again becoming greater with p and, in general, with ^/c. The 
Turbulent mode results are given in figures 17 and 18 at the same 
V/c values as for Cruise, for two of the aircraft. Again the same 
trends are illustrated with higher APF values compared to Cruise, 

SPECTRAL RESPONSE RESULTS 

Power spectral characteristics form the basis for the gust 
response parameters dealt with up to now, such as the A^^ and Ng 
quantities. To gain better understanding of gust response these 
characteristics are now examined in some detail. 

Power spectral density results are summarized in two general- 
ized versions in figures 19 through 34 for the three representa- 
tive flight conditions (Cliit±>/Descent , Low Cruise, High Cruise) 
previously discussed. The basic spectral quantity studied is the 
product of frequency, in Hz, and the acceleration power spectral 
density, in g units squared/Herz . In figures 19 through 26, this 
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parameter is shown, in the first version, normalized by the square 
of rms gust velocity, in ^^/sec^ (ft^/sec^) . The abscissa scale 
in all of the figures is the nondimens ional frequency k divided by 
k for the short period mode frequency as calculated using the pre- 
viously given equation. Thus the scale is ^/ksp» or alternatively 
^/cosp* This course of action was adopted because the peak response 
occurs at the short period natural frequency with A/P off, and 
thus provides a basic reference for comparison of the power spectral 
characteristics of the different airplanes. 

It is seen in figures 19, 21, and 24 for. the controls-fixed 
case that when based on the short period frequencies (which indi- 
vidually do vsry, as noted on the figures, by as much as 87 percent) 
a very close correspondence among peak frequencies is achieved. 

Peak amplitudes differ considerably, and reflect the differences in 
airplane transfer function, f^. Now focusing attention on the 
autopilot effect as seen in a comparison of figures 20, 22 and 25 
with the preceding ones, what is evident is first, the marked gust 
alleviation achieved in all cases and second, the introduction of 
multiple response modes at frequencies both below and above the 
controls-fixed short period frequency. The lower coupled mode is 
fairly flat and wide in character; while the higher one (at the 
A/P on short period frequency of about 2.5 times the controls 
fixed frequency) is sharp and accentuated for two of the four 
airplanes. These high peaks in the spectra can indicate tendencies 
towards resonant-like responses if gust inputs at these frequencies 
are predominant. The major effect of the Turbulent mode (figures 


25 



23 and 26) - at least as seen for two of the four designs - is the 
reduction of sharp peaks when present for the standard A/P settings. 

A similar set of spectral responses is given in figures 27 
through 34 in which the ordinate scale is the acceleration-normalized 
frequency times PSD. This presentation is most convenient for use 
in correlation of flight recorder data since for each spectral 
record obtained a measure of a^n is also obtained. It should be 
noted that there is no difference in curve shape versus frequency 
for a given airplane in this set as compared to the former set. 
Hpwever, the magnitudes of the individual airplane distributions 
relative to each other are altered. Lastly, the results for the 
ANFPSD distributions have been correlated with those obtained from 
flight records for one of the designs. In all essential features 
the comparison proved to be satisfactory. 

CONCLUDING REMARKS 

A simplified approach to the study of aircraft gust response 
allovfs sufficient accuracy for the prediction of vertical c.g. 
acceleration for the frequency bands having the major energy 
content. In commercial jet operations, where autopilot operation 
is in common use, it is necessary to include the autopilot con- 
trol degree of freedom for adequate precision of the gust response 
analysis . 

Autopilot designs follow common trends but individual differ- 
ences in automatic control performance, as seen in a frequency 
analysis, do produce significant differences in the gust response 
characteristics. As used in the four study cases, the pitch-hold 
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mode of the autopilot is of principal importance as compared to 
operation in the altitude-hold and turbulent modes. 

Autopilot effects were such that the gust response parameter 
relating normal acceleration to gust velocity input was reduced 
from 10 to 25 percent depending on the aircraft operating condition. 
The number of zero-level acceleration crossings was increased up 
to 30 percent. Individual aircraft exhibited differences of from 
5 percent to 10 percent, in general, for each of these response 
parameters at the same condition. 

A computer program developed in the study is presented for 
the calculation of gust response characteristics of commercial 
jets with autopilots which are modeled in high order versions. 

Low order simplifications were found to result in significant 
error in the calculation of the response parameters. 
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APPENDIX 


Computer Program 

Donald E. Hawk, Jr 
Villanova University 

The appendix contains a complete presentation of the computer 
program used in the gust response study; the detail given here 
will facilitate reader usage. A particular effort in the develop- 
ment of the program was made to achieve minimum computation time 
for the complex operations involving high-order autopilot simula- 
tion, as well as to afford generality for application to a wide 
variety of autopilot designs. 

The appendix is arranged in terms of the program major ele- 
ments; each beginning with the definitions of all computer names, 
followed by a description of the program element, and ending with 
a flow chart for that portion of the program. The complete program 
is also given. Following is a flow diagram of the overall computer- 
program . 
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Main Program Definitions 


NAPMAX number of highest autopilot mode plus one 

NOftPM omitted autopilot mode 

NFCMAX number of flight conditions 

LALTlto4 lowest flight condition number for a given autopilot 

mode 

NALTlto4 highest flight condition number for a given autopilot 

mode 

C c, wing geometric chord, m (ft) 

CT ctf geometric chord of horizontal tail, m (ft) 


ET 

et, tail moment arm, m (ft) 


KAPPA 

/< , pitch radius of gyration. 

m (ft) 

S 

S, wing area, m2 (ft2) 


ST 

Stf tail area, m2 (ft2) 


AW 

Av 7 , wing aspect ratio 


AT 

At, horizontal tail aspect ratio 

L 

L, turbulence scale, m (ft) 


FC 

fcr cut-off frequency, Hz 


FDLDT 

Fj, control rate effectiveness 

PI 

■ 7T = 3.14159 


KAX(I) 

k, reduced frequency (array) 


HX(I) 

H, altitude, 305m (1000 ft) 

(array) 

XM(I) 

M, Mach number (array) 



RHOX(I) p , air density, kg/m3 (slugs/ft3) (array) 

AOX(I) ag, airplane lift curve slope, rad“l (array) 
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VSX(I) 

ETAX(I) 

FDELX(I) 

EX(I) 

K1 (I) 
K2(I) 
K4(I) 
KY(I) 

AO 

ETA 

PRI 

APRI 

BPRI 

SMLA 

SMLAT 

SM 

MPA 

APM 

WX(1) 

WX (2) 
WX(3) 

IP 

JP 

NFC 

M 

H 


Vg, sonic speed, m/sec (ft/sec) (array) 

= Oe/9a ) , tail downwash derivative (array) 

Fg = Oat/3<S)t/ elevator effectiveness (array) 
e, wing lift anti, m (ft) (array) 

Kl, A/P gain, function of flight condition (array) 

K2, A/p gain, function of flight condition (array) 

K 4 , A/P gain, function of flight condition (array) 

Ky, A/P gain, function of flight condition (array) 

Aq, airplane lift curve slope, rad~l 
= ( 3 e/3 a ) , tail downwash derivative 
prime, parameter of A* and 
A^, wing lift curve slope parameter 
B^, wing lift curve slope parameter 
a, wing lift curve slope, rad“l 
at, horizontal tail lift curve slope, rad“l 
S.M., static margin, percent mean geometric chord 
A/P mode number plus one 
A/P mode 

airplane weight A, 4448N (1000 Ibf) 

airplane weight B, 4448N (1000 Ibf) 

airplane weight C, 4448N (1000 Ibf) 

print-out page number 

causes pages number to be printed 

number of flight condition (used in do- loop) 

M, Mach number of specific flight condition 
H, altitude of specific flight condition 305m (1000 ft) 
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RHO 


VS 

E 

FDEL 

V 

VC 

TVC 

RO 

R1 

COl 

C02 

C03 

GKl 

GK2 

GK4 

GKY 

R2 

SI 

NWT 

W 

MU 

GK3 

GKPHI 


p, air density of specific flight condition, 
kg/m^ (slugs/ft3) 

Vs, sonic speed of specific flight condition, 
m/s (ft/s) 

e, wing lift arm of specific flight condition, ■ 
m (ft) 

Fj, elevator effectiveness for given flight condition 

V, true airspeed, m/s (ft/s) 

V/cf ratio of airspeed to wing chord, sec“l 
2V/c, sec~l 

ro, lift curve slope ratio 

ri, weighted lift curve slope ratio 

e-t/c, tail moment arm to wing chord ratio 

ot/cf tail chord to wing chord ratio 

e/?3, wing lift arm to wing chord ratio 

A/P gain for given flight condition 

A/P gain for given flight condition 

A/P gain for given flight condition 

A/P gain for given flight condition 

r2, horizontal tail arm parameter 

si, horizontal tail distance parameter 

number of the selected weight (used in do-loop) 

W, weight of specific case, 4448N (lOOOlb!') 
n, mass parameter 

A/P gain for given autopilot mode 
A/P gain for given autopilot mode 
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NDIV 


number of divisions to be integrated 


Main Program Description 

The main program performs two tasks - one, the reading of all 
geometric and flight condition dependent data; two, the calculation 
of variables which do not change with k, but do vary with flight 
condition, weight and autopilot mode. 

The first read statement defines the autopilot modes to be 
used, the niomber of flight conditions, and the combination of the 
two for which a case should be run. The statement and an example 
of input are as follows ; 

READ(,5,1)NAPMAX,N0APM,NFCMAX,LALT1,NALT1,LALT2,NALT2, LALT3, NALT3 , 
LALT4,NALT4 
1 FORMAT (1113) 

Col.- 3 6 9 12 15 18 21 24 27 30 33 
53714470025 

This would indicate that: there are 5 autopilot modes APM=0,1,2,3, 

4 (col. 3); that autopilot mode 3 is to be omitted (col. 6); there 
are 7 flight conditions (col. 9); that APM=1 is to be used for 
flight conditions 1 thru 4 (col. 12, 15) , APM=2 for flight conditions 

4 thru 7 (col. 18, 21), and APM=4 for flight conditions 2 thru 5 
(col. 30, 33). 

The next read statement inputs the geometric constants for the 
airplane. The statement and example of input are as follows: 
READ(5,5)C,CT,ET,KAPPA,S,ST,AW,L,FC,FDLDT 

5 FORMAT (life. 0) 
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Col.- 16 12 18 24 30 36 42 48 

24.65 20.10 60.00 38.50 3648. 1338. 7.500 3.790 1000 

54 60 66 

1.200 0.000 

These geometric quantities are then written out by the program in 
a series of write statements. 

The program then calculates and stores the k's (KAX(I)) to be 
used in the detailed print-out. These can be altered by changing 
the do loops and the amount k is indexed. 

The next read statement inputs the specific flight conditions 
parameters such as altitude, Mach number, airplane lift curve 
slope, etc. The following is an example of the statement and 
typical input. 

READ (5,300) HX{1) ,XM(I) ,RHOX(I) ,AOX(I) ,VSX(I) ,ETAX(I) ,FDELX(I) ,EX(I) 
300 FORMAT (8F8. 0) 

Col.- 1 8 16 24 32 40 48 56 64 

10, ,43 .001755 4.96 1078. ,487 .580 4.17 

Next, the gains Kl, K2, K4, Ky which can vary with flight 
condition are read in. The statement and example of input follows; 
READ (5, 310) Kl (I) ,K2 (I) ,K4 (I) ,KY(I) 

310 FORMAT (4F10.0) 

Col.- 1 10 20 30 40 

.379 1.00 1.00 .646 

The quantities that vary with flight condition are then printed - 
that is, the previous two read statement quantities and the static 
margin.. 
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The subroutine AUTO is then called-for each autopilot mode. 
Subroutine AUTO reads all the quantities that vary with autopilot 
mode and writes them out. 

The next read statement is for inputing the weights, WX(1) , 
WX(2),WX(3) in 4448N (1000 Ibf) which corresponds to weights A, 

B, and C respectively. The statement and examples of input follows 
READ(5>1420)WX(1) ,WX(2) ,WX{3) 

420 FORMAT C3F5.0) 

Col.- 1 5 10 15 

260 400. 555. 

These weights are then printed out by a write statement 

All required inputs have been, read into the program at this 
point. Three nested do loops are then used to increment through 
the cases to be run. The outer do loop is for the indexing of 
flight conditions from 1 to NFCMAX. Between this do statement and 
the next are the quantities which vary by flight condition only. 
These quantities are either defined or calculated. The next do 
loop indexes the weight of the aircraft from 1 to 3, or A to C in 
other words. The quantity p' is calculated at this point. The 
inndr-most do loop indexes the autopilot modes. Within this do 
loop the gains that vary with autopilot mode are defined. A series 
of if statements are then used to decide whether this case is to be 
run. This is accomplished by the use of the quantities LALTl thru 
4 and NALTl thru 4 in the if statements. The subroutine INTEG is 
then called (if this case is to be run) and the required quantities 
calculated. These quantities are then written by the subroutine 
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PRINT. The program stops after all do loops have completed in 
dexing the cases. 
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Subroutine AUTO Definitions 


I 

TPT(I, J) 
TPB (I, J) 
TYT(I, J) 
TYB(I, J) 
COND (I) 
CONP (I ) 
CONY (I) 
QUADA(I) 
QUADB(I) 
TDB(I) 
QUAPTA(I) 
QUAPTB (I) 
QUAYTA(I) 
QUAYTB (I) 
QUAPA(I) 
QUAPB(I) 
QUAYA(I) 
QUAYB(I) 
K3 (I) 
KPHI(I) 
NS 
MS 
NPT 
NPB 


autopilot mode 

tau's of phi term, numerator 
tau's of phi term, denominator 
TJy, f tau's of y term, numerator 

tau's of y term, denominator 
C 5 , constant of common term 
, constant of phi term 
Cy, constant of y term 

A, coefficient of common term, numerator 

B, coefficient of common term, numerator 

tau of delta term, denominator 
Ajj^, s 2 coefficient of phi term, nvimerator 
s coefficient of phi term, numerator 
Ay^, s 2 coefficient of y term, numerator 
By^, s coefficient of y term, numerator 
A^ , s^ coefficient of phi term, denominator 
B(j, , s coefficient of phi term, denominator 
Ay, s2 coefficient of y term, denominator 
By, s coefficient of y term, denominator 
A/P gain, function of autopilot mode 
A/P gain, function of autopilot mode 
n, power of s in denominator of phi term 
m, power of s in denominator of y term 
number of tau ' s in numerator of phi term 
number of tau's in denominator of phi term 
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NYT 


m:|inber of tau's fji ?:ivime?:4tor of y term 
nvimber of tau's fn denominator of y term 


NYB 

Subroutine AUTO Description 

The purpose of this subroutine is to read in all quantities 
which are associated with the autopilot control law. All the 
quantities read in, except NPT,NPB,NYT and NYB, can vary with auto- 
pilot mode. The quanitites NPT,NPB,NYT,NYB are fixed at the 
maximum values required. If one autopilot mode has less than the 
maximum number of tau's, zero's are added as needed. The reader 
is referred to the definitions of variables used and the generalized 
control law for an understanding of the following read statements 
and examples of inputs . 

READ (5 , 100 ) NPT, NPB , NYT, NYB 
100 P0RMAT(4I3) 

Col.- 3 6 9 12 
3 3 5 3 

READ (5,110) (TPT(I,J) ,J«1, NPT) 

110 POR^T(10P7.0) 

Col.- 7 14 21 28 35 42 49 56 63 70 

.032 143.0 ,957 

TPB(l,J) ,TyT(I,J) ,andTYB(I,J) read In same manner, 

READ(5,120)COND(I) ,C0NP(I) ,C0NY(I) ,Q0ADA(1) ,QUADB(1) 

120 ?O^T(6F10.Q) 

Col.- 10 20 30 40 50 60 

0.63S 1.00 2.0 3.0 4.0 .135 
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KEAD(5,i23)QUAPTA(I) ,QUAPTB(I) ,QUAYTA{I) ,QUAYTB(I) 

123 FORMAT (4F1Q,0) 

Col,- 10 20 30 40 

i00117 0,5263 0,0 0.0 

QUAPA(I) ,QUAPB(1) ,QUAYA(I) ,QUAYB(I) read in a similar manner. 
READ (5, 130) K3 (I) ,KPHI(I) 

130 FORMAT (2F1 0.0) 

Col.- 10 20 

1,0 1,0 

READ (5, 135) NS (I) , MS (I) 

135 FORMAT (213) , 

Col.- 36. 

1 0 

The quantities read in are then written out with the proper des 
criptions . 
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KC 

DELTA 

ANS 

ANSI 

EOT 

TOP 

BIGK 

KPRIME 

KO 

NO 

GLEG15 
FUN 
FUNl 



Subroutine INTEG Definitions 

kc, reduced cut-off frequency 

width of interval to be integrated 

summation of values GLEG15 

summation of values GLEG15 

lower limit of integration interval 

upper limit of integration interval 

K, gust response parameter 

K , gust response parameter 

kQ, ratio of gust response parameters 

Nq, number of positive-slope crossings per second 

of one-g load level due to gusts 

function which equals value of integral 

names of subroutines which generate required 

values for integration 

Subroutine INTEG Description 


The purpose of this subroutine is to perform the integrations 
required for K and K^, The upper limit of the integral, kg, is 
calculated and the interval from 0 to the upper limit, kc » is 
divided into 10 equal divisions. Each division is then integrated 
twice, once for the K calculation and once for the calculation. 
This integration is performed using a 15 point Gauss-Legendre Quad- 
rature provided by the statement function GLEG15 and the integrals 
of the divisions are summed to provide the total integral from 0 
to k^. Then the quantities K, K^, kQ, and Nq are calculated using 
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the previously found integrals. 










GLEG15 Definitions, evaluates integrals G(x)dx 

W, weighting for results of function integrated, G(x) 
R, weighting for location of point in interval 

A, lower limit of integral 

B, upper limit of integral 

C, half the distance from A to B 

D, average of A and B 
k, reduced frequency (array) 
value of integral in interval A to B 
external subroutine to calculate the value of 
G(x) at specified points, FUN, FUNl are the 
external subroutines. 

Function GLEG15 Description 

The function GLEG15 evaluates the integral of a function from 
the lower limit. A, to the upper limit, B. The name of the func- 
tion, and limits, A and B, appear in the function statement in 
subroutine INTEG. The two external subroutines which are used in 
conjunction with GLEG15 are FUN AND FUNl. Subroutine FUN calculates 
the function values for the integral used in the calculation of K. 
Subroutine FUNl calculates the values for the integral used in the 
calculation of K^, After the integral is calculated over the 
interval using a 15 point Gauss-Legendre Quadrature the value is 
returned to subroutine INTEG, 


Function 

W(l) to W(8) 
R(l) to R(8) 

A 

B 

C 

D 

KAG(I) 

GLEG15 

FUNCT 
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J 

J 
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Subroutine PUN and FUNl Definitions 


KAG(J) 

KAY 

FI 

TLC 

TLCK 

PHISIG 

F(J) 

CALFl(J) 


k, reduced frequency (array) 

k, reduced frequency value 

fl, airplane generalized transfer function 

2L/c/ twice the ratio of turbulence scale to 

wing chord 

parameter in 
2 

'^’w/ai , ratio of vertical gust velocity PSD to rms 
vertical gust velocity square 
values of points to be used in Gauss Legendre 
Quandrature, GLEG15 (array) 

fl, stored value of fj calculated in subroutine FUN 
to be used again in subroutine FUNl (array) 


Subroutine FUN Description 

This subroutine is used to obtain the function values necessary 
in the calculation of the integral for K. The calculations are 
performed at the 15 values of k specified by the Function GLEG15. 
Subroutine CALC is called upon to obtain the value of the airplane 
generalized transfer function fi at each specified k. These values 
are stored for use in subroutine FUNl. The function values to be 
used in the integration are calculated and returned to GLEG15 as 
the quantity F(15). 


50 



Subroutine FUNl Description 

This Subroutine is used to obtain the function ^ values nec- 
essary in the calculation of the integral for k'". Calculations 
are performed at the 15 values of k specified by the function 
GLEG15. The basic quantity to be calculated is the airplane 
generalized transfer function f^. Since FUN had already obtained 
these values at the same k's as FUNl required them, the values of 
fl are stored in FUN and used again in FUNl. The function values 
to be used in the integration are calculated and returned to 

I 

GLEG15 as the quantity F(15), 


[ 

1 
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Subroutine CALC Definitions 


KAY 

OMEGA 

SOP 

DN 

DD 

PN 

PD 

YN 

YD 

STAB 

TAIL 

ABPHI 

ABY 

PARI 

PAIl 

PBRl 

PBIl 

Cl 

CIBAR 

PAR2 

PAI2 

PBR2 

PBI2 

C2 

C2BAR 


k, reduced frequency specific value 
o, (2V/C)k, response frequency 
s, Laplace operator 

denominator of common term of A/P equation 

numerator of common term of A/P equation 

numerator of phi term of A/P equation 

denominator of phi term of A/P equation 

numerator of y term of A/P equation 

denominator of y term of A/P equation 

imaginary term of TAIL 

coefficient of ABPHI and ABY 

airplane phi coefficient parameter 

airplane y coefficient parameter 

component of a^, real 

component of a^, imaginary 

component of bj, real 

component of b^, imaginary 

cx airplane coefficient, real part 

Cl airplane coefficient, imaginary part 

component of a 2 , real 

component of a 2 , imaginary 

component of b2, real 

component of b2, imaginary 

C 2 airplane coefficient, real part 

C 2 airplane coefficient, imaginary part 



PAl 

P.B1 

PA2 

PB2 

AlCOM 

A2COM 

BICOM 

B2COM 

CICOM 

C2COM 

YTOP 

YBOT 

Y 

YDOT 

YDTC 

YDOTSQ 

BETASQ 

F2 

F3 

FI 


component of a^, complex 

component of b^, complex 

component of a 2 » complex 

component of b 2 f complex 

aj + a 2 airplane coefficient, complex 

32 + 32 airplane coefficient, complex 

bi + b^ airplane coefficient, complex 

b 2 + b 2 airplane coefficient, complex 

Cj^ + Cq^ airplane coefficient, complex 

C 2 + C 2 airplane coefficient, complex 

niomerator of vertical displacement equation 

denominator of vertical displacement equation 

y, vertical displacement, m (ft) 

• 

y, vertical velocity real component, m/s (ft/s) 

• 

yi, vertical velocity imaginary component, m/s (ft/s) 

|y|2, vertical velocity's magnitude squared, m^/s^ (ft^/s^ 
Mach number correction, Glauret factor square 
f 2 , airplane transfer function correction factor 
f^, airplane transfer function correction factor 
fq, airplane generalized transfer function 


Subroutine CALC Description 


The purpose of this subroutine is to obtain the value of the 
airplane generalized transfer function fq, for a specified reduced 
frequency k. This required the calculation of the quantity lyp, 
which is the absolute magnitude of the vertical velocity vector 
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square. This, in turn, would have necessitated determining the 
vertical displacement y from a 5x5 matrix. The matrix is com- 
posed of the two equations of motion, the two downwash equations, 
and the autopilot control law. The solution of this set of simul- 
taneous equations for y was found in closed f oann to afford a three- 
fold reduction in computer time. The solution was obtained by 
reducing the 5x5 matrix to a 2x2 matrix through successive substi- 
tutions and then applying Crammer ’ s rule . The subroutine obtains 
the coefficients of the 2x2 matrix, which are termed a^^, b^^, C]^, 
a 2 , b 2 , and C 2 « These coefficients are complex and the program 
recognizes them as such. The major part of the subroutine involves 
the calculation of components of these coefficients, both real and 
imaginary. Once y is obtained the calculation of |y|2 is performed. 
The airplane transfer function fi is then calculated and returned 
to the calling subroutine. 
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Subroutine CALC Flow Diagram 


Do J- 


Do 


Do J' 


Do 




























Equations for Subroutine CALC 



ioj 


PN = C^K^ s^+b(j,^ s+ 1 ) (X(j,^s+ 1 ) 2 ^. . . s+l)^PY 

PD = sn(a^^s2+b^^s+l) (t^^s+1)i... 

YN = CyKy (ay^ s2+by^ s+1) (xy^s+1) 2 ^. . . (xy^s+1)^^^ 

YD = (ay^s2+by^S + l) (Xy^S + 1) . . (Xy^S+1)^^^ 

DN = s+1 

K4 

DD = C K 1 K 2 K 3 (as2+bs+l) 

STAB = F 5 ft k 
c 

TAIL = -Fj -{F 5 ft k)i = -Fg - STABi 

c 

ABPHI = TAIL (5£) (^) 

DN PD 

ABY = TAIL (25) (M) 

DN YD 

PARI = -2ia(2V)k2 
c 

PAIl = (IZ) [ro+ri(l-T^)']k 

PBRl = -[ro+ri(l-r^^ 

PBIl ^ [ro(|^ - D-ri r 2 ]k 

Cl = d = ro + ri [cos kS 2 ^ - r\"] 

CIBAR = Cl = -ri(sin kSi)i 
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PAR2 


0 


PAI2 = (|Z) ro - ^ ri (l-r^ijk 
PBR2 = - r® ro - ft xi (1-T^) + 4uk^k^l 

Lc c 

PBI2 = r® rg (2£ _ 1) + ^ k 

*-C C c 

C2 = C 2 = ® rg - ft (cos kSn 
c c 

C2BAR = C 2 = ft (sin kSj) 
c 

PAl = -2Mak2 + (^) [ro + ri (1 -t^)“] ki 
c c , 

= PARI + PAIli 
PBl = -[rg + 

= PBRl + PBIli 

PA2 = 0 + (^) p. rg - ^ ki 

= PAR2 + PAl2i 

2 

PB2 = -[| rp - ft ri (l-ri) + 4-w ^2 ^2j 


[I -^0(21- 


1) + ®t rir^ ki 


= PBR2 + PBI2i 


AIGOM = ai + ai = riAy - 2y(^)k2 + ^riBy + (^)[^ro + ri(l-r^)^ k^i 

= r^CAy + Byi) - 2y(^)k2 + (If ) |j:g + ri(l-T\)]j ki 


= RI ABY + PAl 
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A2C0M 


BICOM 


B2COM 


CICOM 


C2COM 


YTOP = 


= a 2 + a 2 


_e 

c 


t riAy + riBy + (2V) fe ro - ft ri i 

V ^ O o o 

ft ri(Ay+Byi) + (|V) ro - ft rid-r;^ ki 


C -G C 

-COl Rl ABY + PA2, (COl=et/c) 


= + bi = 


tiA^ “ C^O + ri(l-T^j] + “l^riB,,, + [ro(|^ - l)-rir2']kji 
= ri(A^+B<{,i) - [^ro + (1-r^^ + [jo 


= Rl ABPHI + PBl 


b2 + b2 = -ft riA* - ro - ft ri(l-i^) + 4y k^l 
c Lc c c J 

+ l^-ftiXiB^ + j”® ro (^ - 1) + ft rir^kj- i 

= -ft ri(A^+B(|)i) - [~£ ro - ft ri(l-n) + 4 m-^2 
c tc q . ® -* 

= -COl Rl ABPHI + PB2 

cj + = ro + ri (cos kS^-n)- r^ (sin kS^)! 

= Cl + ClBARi 


0 0. 

= C 2 + C 2 = — ro - _t r]^(cos kSj^ - r\) + r^ (sin kSi)i 


= C2 + C2BARi 


(ci+ci) (b2+b2) - (C2+C2) (bi+bi) 
CICOM B2COM - C2COM BlCOM 
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YBOT = (ai+aj) (b2+b2) - (a2+a2) (bi+bi) 


= AlCOM B2COM - A2COM BlCOM 

Y = y (ci+ci) (b 2 +b 2 ) - (C2+C2) (bi+bj ) 
(ai+ai) (b 2 +^ 2 ) “ (a 2 +a 2 ) (bx+Fi) 


= YTOP 
YBOT 

YDOT = y = icjy = i k)y = (y)r + (y) i 

c 

• • • 

YDTC = conjugate of y = (y)r - (y)i 

ly |2 = (y) - (y)^ 2 

= [(y)r + (y)i^[(y)r - ( y ) i \ 

= (YDOT) (YDTC) 

BETASQ = = 1-m2 


F2 = f2 = 1_ 

1 + 1.5-rrk + ir2 Mk2 

3^ 3^ 


F3 = fg = 1 

1 + .55 Awk 

FI = fi = 4^x2 y.2 |y|2 f 2 f 3 
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Subroutine PRINT Definitions 


IPRT 


KAY 
KAX(J) 
FIX (J) 
TLC 

TLCK 
YDTSQ ( J) 
F2X(J) 
F3X (J) 
FISIG(J) 

JP 

KP 

IP 

ASIG 

ASIGVC 

F2F3 

FUNFl 

FUNIFI 

GFPSD 


controls whether detailed print-out is to be 
obtained for that case 

The following arrays are found for detailed print-out 
at the selected reduced frequencies. 

k, selected reduced frequency 
k, reduced frequency (array) 

fl, airplane generalized transfer function (array) 

2L/c, twice the ratio of turbulence scale to wing 
chord 

2 

parameter in ‘f’w/ai 

|yp, absolute value of velocity square (array) 

f2/ airplane transfer function correction factor (array) 

f 3 , airplane transfer function correction factor (array) 
2 

<{iw/a]_ , ratio of vertical gust velocity PSD to rms 

vertical gust velocity square (array) 

causes page number to be printed 

causes page nximber to be printed 

page number of print out 

Aa, gust response parameter 

Aa/(V/fc), generalized gust response parameter 
f2 f 3 , airplane transfer function correction factor 
values generated to find K 
values generated to find 

GFPSD, generalized frequency times power spectral 
density 
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GNFPSD 


GNFPSD, gust normalized frequency times power 
spectral density 
ANFPSD ANFPSD, acceleration normalized frequency times 

power spectral density 

Subroutine PRINT Description 

The purpose of this subroutine is to print out the results 
for the cases run. The standard printout is performed for all 
cases. The detailed printout is done only for the medium weight, 

B, cases. The standard printout is of the following quantities: 
case number, M, h, V, V/c, y , kc, aQ, a, a^, rg, r^, r 2 , s^, K, 

K^, Ao, Aa/(V/c)', kg, Nq . The detailed printout are the values of 
the quantities: \y\^r ^1' r fi(‘(’w/ax )/ k2 f^^ ‘f’w/ai^, 

GFPSD, GNFPSD, ANFPSD calculated for specific values of k. The 
values of k are specified in the main program and run for three 
decades, 10“^ to 1, The purpose of the detailed printout is to 
generate the PSD graphs. Subroutine CODE, which is used in con- 
junction with subroutine PRINT, prints out the case number when 
called upon. Once the required printout for the case is performed 
the program returns to the Main Program. 
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Subroutine CODE Flow Diagram 




Program Listing 


* $$ JOB JNH-HAHK,CLASS-A 

* $$ PRT CLASS-H 

// JOB HAWK, 3610051010 
// EXEC PROC-FRTCLG 

IMPLICIT REAL*8(A-H,K,0-Z) 

REAL*8 MO,M,NRSQ,NISQ,L 
INTEGER APM 

COMMON/ INP/H , RHO , VS , AO , ETA , FDEL , FDLDT , W, C , CT , ET , S , ST , AW , AT , FC 

I, KAX(60) , KAPPA, RAY, M,L 
COMMON/FLT/APM, NFC , NWT 

COMMON/CALCU/ V, VC , PI , APRI , BPRI , SMLA , SMLAT , R0 , R1 , E , COl , C02 , C03 , R2 , S 

II, MO 

COMMON/ INT/BIGK , AS IG , FI , KPRIME , KC , KO, NO 
COMMON/PRT/YDOTSQ , F2 , P3 , IP , JP 
COMMON/FAST/TVC,GKl ,GK2 ,GK3 ,GK4 ,GKPH1,GKY 
COMMON/DIV/NDIV 

COMMON/GEN/TPT(4,10) ,TPB(4,10) ,TYT(4,10) ,TYB(4,10) ,COND(4) ,CONP(4) 
l,CONY(4) ,QUADA(4) ,QUADB(4) ,TDB(4) ,QOAPA(4) ,QUAPB(4) ,QUAYA(4) ,QUAYB 
2(4) ,QOAPTA(4) ,QUAPTB(4) ,QUAYTA(4) ,QUAYTB(4) ,NS(4) ,MS(4) 
COMMON/GEN2/K3(4) ,KPHI(4) ,NPT,NPB,NYT,NYB 

DIMENSION XM(20) ,HX(20) ,RBOX(20) ,VSX(20) ,AOX(20) ,ETAX(20) ,FDELX(20 
1) ,EX(20) ,K1(20) ,K2(20) ,K4(20) ,KY(20) ,HX(3) 

READ (5,1) NAPMAX , NOAPM, NFCMAX , LALTl , NALTl , LALT2 , NALT2 , LALT3 , NALT3 , L 
1ALT4,NALT4 

READ(5,5)C,CT,ET,KAPPA,S,ST,AW,AT,L,FC,FDLDT 

PI-3.14159 

WRITE (6, 10) 

WRITE(6,15)C 
WRITE(6,20)CT 
WRITE (6, 25) KAPPA 



WRITE (6 ,30) ET 
WRITE(6,35)S 
WRITE(6 ,40) ST 
WRITE(6,45)AW 
WRITE (6, 50) AT 
WRITE (6, 55) L 
WRITE(6,60)FC 
WRITE(6,65)FDLDT 
KAX(1)=0.001 
DO 100 1 = 2,19 
KAX(I)=KAX(I-1)+.001 
100 CONTINUE 

DO 120 1-11,28 
KAX(I)-KAX(I-l)+.005 
120 CONTINUE 

DO 140 1=29,33 
KAX(I)=KAX(I-l)+.02 
140 CONTINUE 

DO 160 1-34,41 
KAX(I)=KAX(I-l)+.l 
160 CONTINUE 

WRITE (6, 250) 

DO 390 I-1,NFCMAX 

READ(5,300)HX(I) ,XM(I) ,RHOX(I) ,AOX(I) ,VSX(I) ,ETAX(I) ,FDELX(I) ,EX(I 

390 CONTINUE 

DO 400 I-1,NFCMAX 
READ(5,310)K1(I) ,K2(I) ,K4(I) ,KY(I) 

AO-AOX(I) 

ETA-ETAX(I) 

PRI-.9*PI/((1./AT)-(1./AW)) 

APRI-AO- (1 .+ ( ST/S) * ( 1 .-ETA) ) *PRI 
BPRI-AO*PRI 



SMLA«(APRl/2.) +DSQRT( (APRl/2.) **2+BPRI) 

SMLAT-SMLA/ (1 .+SMLA/PRI) 

SM- (SMLA/C/AO) *100 . * (SMLAT*ST* (1 ,-ETA) *ET/SMLA/S-EX ( I) ) 

WRITE(6,340) I,XM(I) ,HX(I) ,RHOX(I) ,VSX(I) ,AOX(I) ,ETAX(I) ,FDELX(I) ,E 
1X(I) ,SH,K1{I) ,K2(I) ,K4{I) ,Ky(I) 

400 CONTINUE 

DO 410 MPA-1, NAPMAX 
APH-HPA-1 

IP (APM. EQ. NOAPM) G0T041 0 
CALL AUTO 
410 CONTINUE 

READ(5,420)WX(1) ,WX(2) ,WX(3) 

WRITE(6,S00) 

WRITE(6,550)WX(1) ,WX(2) ,WX{3) 

IP-0 

JP-0 

DO 600 NFC-1, NFCMAX 
M-XH(NFC) 

H-HX(NFC) 

RHO-RHOX(NFC) 

VS-VSX(NFC) 

AO-AOX(NFC) 

ETA- ET AX (NFC) 

E- EX (NFC) 

FDEL-FDELX(NFC) 

V-M*VS 

VC-V/C 

TVC-2.*VC 

PRI-.9*PI/((1./AT)-(1./AW)) 

APRI-AO-(l.+(ST/S) *(1.-ETA) ) *PRI 
BPRI-AO*PRI 

SMLA- (APRI/2. ) +DSQRT ( (APRI/2. ) **2+BPRI) 

SMLAT- SMLA/ ( 1 . +SMLA/PRI ) 



R0=SMLA/AO 

R1 » ( SMLAT/AO ) * ( ST/S) 

C01=ET/C 

C02=CT/C 

C03«=E/C 

GK1=K1(NFC) 

GK2=K2(NFC) 

GK4-K4(NFC) 

GKY-KY(NFC) 

R2 = ( 2 . *C01 ) +C02+ ( ETA* ( 2 . *C03 -1 . ) ) 

S1«2.*C03+2.*C01-1.+C02 

DO 650 NWT-1,3 

W=WX(NWT) 

MU»2.*W*1000./(32.2*AO*RHO*S*C) 

DO 700 MPA=1,NAPMAX 
APM=MPA-1 

IF(APM.EQ.NOAPM)GOTO 700 

GK3=°0.0 

GKPHI-0.0 

IF (APHi NE. 0) GK3=K3 ( APM) 

IF (APM. NE. 0) GKPHIfKPHI (APM) 

GKY«0. 

IF (APM. EQ. 2) GKY=KY (NFC) 
IF(NFC.LT,LALTl.AND.APM.EQ.l) GO TO 700 
IF(NFC.GT.NALTl.AND.APM.EQ.l) GO TO 700 
IF(NFC.LT.LALT2.AND.APM.EQ.2)GO TO 700 
IF(NFC.GT,NALT2.AND.APM.EQ.2)GO TO 700 
IF(NFC.LT.LALT3.AND.APM.EQ.3)GO TO 700 
IF(NFC.GT.NALT3.AND.APM.EQ.3)GO TO 700 
IP(NFC.LT.LALT4.AND.APH.EQ.4)GO TO 700 
IF(NFC.GT.NALT4.AND.APM.EQ.4)GO TO 700 
800 CONTINUE 


NDIV=10 



CALL IMTEG 


CALL PRINT 
700 CONTINUE 
650 CONTINUE 
600 CONTINUE 
STOP 

1 FORMAT (1113) 

5 FORMAT (11F6.0) 

10 FORMAT (IX, '*****************************•,/, lx, 'AIRPLANE GEOMETRIC 
1 CONSTANTS' ,/, lx, '*****************************' ,/) 

15 FORMAT (IX, 'WING GEOMETRIC CHORD,C* ,6X,P5.2, ' FT.') 

20 FORMATdX, 'TAIL GEOMETRIC CHORD, CT' ,5X,F5. 2, ' FT.') 

25 FORMATdX, 'PITCH RAD. OF GYRATION, KAPPA ',F4.1,' FT.') 

30 FORMATdX, 'TAIL MOMENT ARM,ET* ,10X,F5.1, ' FT.') 

35 FORMATdX, 'WING AREA, S' ,17X,F5.0 , ' SQ.FT.') 

40 FORMATdX, 'TAIL AREA, ST' ,16X,F5.0,' SQ.FT.') 

45 FORMATdX, 'WING ASPECT RATIO,AW' ,9X,F4.2) 

50 FORMATdX, 'TAIL ASPECT RATIO, AT' ,9X,F4. 2) 

55 FORMATdX, 'TURBULENCE SCALE, L',10X,F5.0,' FT.') 

60 FORMATdX, 'CUT-OFF FREQUENCY, FC ,9X,F4. 2,' HZ. ') 

65 FORMATdX, 'ELEVATOR RATE EFFECTIVENESS •,F3.1,///) 

250 FORMATdX, 'NFC ,3X, 'M' ,4X, 'H' ,5X, 'RHO* ,5X, *VS' ,5X, 'AO' ,3X, 'ETA' ,2X 
1,'FDEL' ,2X, ' E ',3X,'S.M. ',4X,'GK1',5X,'GK2',5X,'GK4',5X,'GKY') 

300 FORMAT(8F8.0) 

310 FORMAT (4F10.0) 

340 FORMAT(1X,I2,2X,F4.2,2X,F3.0,1X,F8.6,2X,F5.0,2X,F4.2,1X,F5.3,1X,F5 
1.3,1X,F4.2,2X,F4.1,4(2X,F6.3) ,/) 

420 FORMAT(3F5.0) 

500 FORMAT(///, IX, 'GROSS WEIGHT',/) 

550 FORMAT (4X, ' A- ' ,F5.0/,4X,'B-' ,F5.0/ ,4X, 'C-' ,F5.0//) 

END 

SUBROUTINE AUTO 
IMPLICIT REAL*8(A-H,K,0-Z) 



REAL*8 MU,M,NRSQ,NISQ,L 
INTEGER ARM 

COMMON/GEN/TPT(4,10) ,TPB(4,10) ,TYT(4,10) ,TYB(4,10) ,C0ND(4) ,C0NY(4) 
1,C0NP(4) ,QUADA(4) ,QUADB(4) ,TDB(4) ,QUAPA(4) ,QUAPB(4) ,QUAYA(4) , QUAYS 
2(4) ,Q0APTA(4) ,QUAPTB(4) ,QUAYTA(4) ,Q0AYTB(4) ,NS(4) ,HS(4) 
COMHON/GEN2/K3(4) ,KPHI(4) ,NPT, NPB, NYT, NYB 

COMMON/ INP/H , RHO , VS , AO , ETA , FUEL , FDLDT , W, C , CT , ET , S , ST , AW, AT , PC 
1,KAX(60) , KAPPA, KAY, M,L 
COMMON/FLT/APM, NFC , NWT 
IF(APM.EQ.0)WRITE(6,460) 

IF (APM.EQ.0) RETURN 

IF ( APM. EQ . 1 ) WRITE (6,461) 

IF(APM.EQ.2)WRITE(6,462) 

IF (APM. EQ. 3) WRITE (6, 463) 

IF (APM. EQ. 4) WRITE (6, 464) 

I=APM 

DO 105 J-1,10 
TPT(I,J)-0.0 
TPB(I,J)-0.0 
TYT(I,J)=0.0 
TYB(I,J)«0.0 
105 CONTINUE 

READ ( 5 , 1 0 0 ) NPT , NPB , NYT , NYB 
READ(5,110) (TPT(I,J) ,J»1,NPT) 

READ(5,110) (TPB(T,J) ,J-1,NPB) 

READ(5,110) (TYT(I,J) ,J«1,NYT) 

READ(5,110) (TYB(I,J) ,J=1,NYB) 

READ(5,120)COND(I) ,CONP(I) ,COHY(I) ,QUADA(I) ,QUADB(I) ,TDB(I) 
READ(5,123)QUAPTA(I) ,QUAPTB(I) ,QUAYTA(I) ,QUAYTB(I) 
READ(5,125)QUAPA(I) ,QUAPB(I) ,QOAYA(I) ,QUAYB(I) 

READ(5,130)K3(I) ,KPHI(I) 

READ(5,135)NS(I) ,MS(I) 

WRITE(6,140) (TPT(I,J) ,J-1,NPT) 
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WRITE(6,145) (TPB(I,J) ,J-1,NPB) 

WRITE (6“,15B) (TYT(I,J) ,J»1,NYT) 

WRITE(6,155) (TYB(I,J) ,J-1,NYB) 

WRITE(6,160)TDB(I) 

WRITE(6,161)QUADA(I) ,QUADB(I) 

WRITE(6,162)QUAPA(I) ,QUAPB(I) 

WRITE(6,172)QUAPTA(I) ,QOAPTB(I) 

WRITE(6,163)QUAYA(I) ,QUAYB(I) 

WRITE (6,173 )QUAYTA( I) ,QUAYTB(I) 

WRITE(6 ,165)COND(I) ,CONP(I) ,CONY(I) 

WRITE(6,175)K3(I) ,KPHI(I) 

WRITE(6 ,170)NS(I) ,MS(I) 

RETURN 

100 FORMAT(4I3) 

110 FORMAT (10F7.0) 

120 FORMAT (6F10.0) 

123 FORMAT(4F10.0) 

125 FORMAT (4F10.0) 

130 FORMAT (2F10.0) 

135 F0RMAT(2I3) 

140 FORMAT ( IX, /, 9X, 'TAUS FOR NUMERATOR OF VPHI ' ,4X ,10F8 . 4) 

145 FORMAT(lX,/,9X, 'TAUS FOR , DENOMINATOR OF VPHI ' , 2X , 10F8 . 4) 

150 F0RMAT(1X,/,9X, 'TAUS FOR NUMERATOR OP V Y ',4X,10F8.4) 

155 FORMAT{lX,/,9X, 'TAUS FOR DENOMINATOR OF V Y ',2X,10F8.4) 

160 FORMAT ( IX, /, 9X, 'TAU FOR DENOMINATOR OF COMMON TERM',F10.6) 

161 FORMAT ( lx, /,9X, 'COMMON QUADRATIC COEFFICIENTS : A- • ,F10 . 6 ,2X ,' B= ' ,P 

110 . 6 ,/) 

162 FORMAT (9X, 'DENOMINATOR PHI QUADRATIC COEFFICIENTS: A- ' ,F10. 6 ,2X , ' B 
1*' ,F10.6,/) 

163 FORMAT (9X, 'DENOMINATOR Y QUADRATIC COEFFICIENTS: A- ' ,P10 . 6 ,2X, ' B» 
l',F10.6,/) 

165 FORMAT (9X, 'FIXED GAIN OF COMMON TERM- ' ,F10 .6 ,// ,9X, ' FIXED GAIN OF 
IPHI TERM-' ,5X,F10. 6,//, 9X, 'FIXED GAIN OF Y TERM- ' ,5X ,F10 .6 ,/) 



172 FORMAT (9X, 'NUMERATOR PHI QUADRATIC COEFFICIENTS: A= ' ,F10 . 6 ,2X , ' B= ' 
1,F10.6,/) 

173 FORMAT (9X, 'NUMERATOR Y QUADRATIC COEFFICIENTS: A« ' ,F10 . 6 ,2X , ' B= ' , 
1F10.6,/) 

175 FORMAT (9X, 'VARIABLE GAIN GK3-' ,F10. 6 ,// ,9X ,' VARIABLE GAIN GKPHI-', 
1F10.6,/) 

170 FORMAT (9X, 'EXPONENTS: N-',I2,' M-',I2,/) 

460 FORMAT (/////, 2X, 'APM ' , ' 0 ' ,2X, ' A/P OFF',/) 

461 FORMAT (2X, 'APM ', '1 ' ,2X, 'PITCH HOLD',/) 

462 FORMAT (2X, 'APM ' , '2' ,2X, 'PITCH+ALTITUDE HOLD',/) 

463 FORMAT (2X, 'APM ',' 3 ', 2X , 'TURBULENT MODE',/) 

464 FORMAT (2X, 'APM ' , ' 4 ' , 2X , ' PITCH+RATE HOLD',/) 

END 

SUBROUTINE INTEG 
IMPLICIT REAL*8(A-H,K,0-Z) 

REAL*8 NO,M 
INTEGER APM 

COMMON/INP/H,RHO,VS,AO,ETA,FDEL,FDLDT,W,C,CT,ET,S,ST,AM,AT,FC 

I, KAX(60) , KAPPA, KAY, M,L 
COMMON/FLT/APM, NFC , NWT 

COMMON/CALCU/V, VCjPI ,APRI ,BPRI,SMLA, SMLAT, R0 ,R1 ,E,C01 ,C02 ,C03 ,R2 ,S 

II, MU 

COMMON/ INT/BIGK , AS IG , FI , KPRIME , KC , KO > NO 

COMMON/PRT/YDOTSQ, F2 ,F3 , IP, JP 

COMMON/DIV/NDIV 

EXTERNAL FUN,FUN1 

KC»FC*PI/VC 

DELTA- KC/NDIV 

ANS-0.0D0 

ANSl-0.000 

BOT-0.0D0 

TOP-DELTA 

DO 10 I-1,NDIV 
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ANS=AHS+GLEG15 (BOT,TOP,FUN) 

ANS1-ANS1+GLEG15 (BOT, TOP, FUND 
TOP-TOP+DELTA 
BOT-BOT+DELTA 
10 CONTINUE 

BIGK«DSQRI(ANS) 

KPRIME-DSQRI(ANSl) 

KO-KPRIME/BIGK 

NO-VC*KO/PI 

RETURN 

END 

FUNCTION GLEG15(A,B,FUNCT) 

IMPLICIT REALMS (A-H,K,0-Z) 

C 

C FUNCTION GLEGIS EVALUATES INTEGRALS OF THE FORM 
C 

C B 

C INTEGRAL G(X) OX 

C A 

C : 

C WHERE A AND B ARE FINITE USING IS POINT GADSS-LEGENDRE QUADRATURE. 

C 

C USE 
C 

C VALUE-GLEG1S(A,B,FUNCT) 

C 

C A THE LOWER LIMIT OF INTEGRATION. MUST BE FINITE. 

C 

C B THE UPPER LIMIT OF INTEGRATION. MUST BE FINITE. 

C 

C FUNCT..THE NAME OF AN EXTERNAL SUBROUTINE TO CALCULATE THE 

C VALUE OF G(X) AT SPECIFIED POINTS. 

C SUBROUTINE FUNCT(F,N) 
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C F..ON INPUT THE N VALUES OF X 

C ON OUTPUT THE N VALUES G(X) 

C F(I)I 1=1,. ..,N 

C N..THE NUMBER OF FUNCTION EVALUATIONS 

C THE NAME GIVEN FUNCT MUST APPEAR IN AN EXTERNAL 

C STATEMENT IN THE CALLING PROGRAM 

C 

C GLEG15 IS AVAILABLE IN FORTRAN 
C 

C RELEASED 25 AUGUST 1975 
C 

C PROGRAMMED BY S. BAUDENDISTEL COMPUTER SCIENCES CORPORATION 
C 

COMMON/FUNY/F(15) ,KAG(15) ,CALF1(15) 

DIMENSION W(8) ,R(8) 

W(l)=. 202578241925561 
W(2)=. 198431485327111 
W(3)=. 186161000015562 
W(4)=. 166269205816993 
W(5) -.139570677926154 
W(6) -.107159220467171 
W(7)-. 07 0366047488108 
W(8)-. 030753241996117 
R(1)-0.0D0 

R(2)=. 201194093997434 
R(3)-. 394151347077563 
R(4)-. 570972172608538 
R(5) =.72441773136017 
R(6) =.848206583410427 
R(7) -.937273392400705 
R(8)=. 987992518020485 
C=(B-A)*.5 
D=(B+A)*.5 
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KAG(l) >0 
DO 10 1-2,8 
S-R(I)*C 
KAGCD-O+S 
J-I+7 

KAG(J) -D-S 
10 CONTINUE 
CALL FUNCT 
D-W(l) *F(1) 

DO 20 1-2,8 
J-I+7 

20 D-D+W(I)*(F(I)+F(J)) 

GLEG15-C*D 

RETURN 

END 

SUBROUTINE FUN 

IMPLICIT REAL*8(A-H,K,0-Z) 

REAL+8 M 
REAL* 8 L 
INTEGER APM 

COMMON/ INP/H , RHO , VS , AO , ETA , FDEL , FDLDT , W, C , CT , ET , S , ST , AH, AT , FC 

I, KAX(60) , KAPPA, KAY, M,L 
COMMON/FLT/APM , NFC , NWT 

COMMON/ CALCU/V, VC , PI , APRI , BPRI , SMLA , SMLAT , R0 , R1 , E , COl ,C02 ,C03 , R2 , S 

II, MU 

COMMON/ INT/BIGK , AS IG, FI ,KPRIME , KC , KO , NO 
COMMON/PRT/YDOTSQ , F2 , F3 , IP , JP 
C0MM0N/FUNY/F(15) ,KAG(15) ,CALF1(15) 

DO 10 J-1,15 
. KAY-KAG(J) 

CALL CALC 

CALF1(J)-F1 

TLC-2.+L/C 


79 



TLCK-(1.339*TLC*KAY)**2 

PHISIG= (TLC** ( 5 ,/3 . ) ) * (1 . + ( 8./3 . ) *TLCK) / ( (1 .+TLCK) ** (11 ./ 6 .)) 
F(J)-F1*PHISIG 
10 CONTINUE 
RETURN 
END 

SUBROUTINE FUNl 
IMPLICIT REAL*8(A-H,K,0-Z) 

REAL*8 M 
REAL* 8 L 
INTEGER APM 

COMMON/ INP/H , RHO , VS , AO , ETA , FDEL , FDLDT , W , C , CT , ET , S , ST , AW , AT , FC 

I, KAX(60) , KAPPA, KAY, N,L 
COMMON/FLT/APM, NFC , NWT 

COMMON/ CALCU/V , VC , PI , APRI , BPRI , SMLA , SMLAT , R0 , R1 , E , COl , C02 , C03 , R2 , S 

II, MU 

COMMON/ INT/BIGK , AS IG , FI , KPRIME , KC , KO , NO 
COMMON/ PRT/YDOTSQ , F2 , F3 , IP , JP 
COMMON/FUNY/F(15) ,KAG(15) ,CALF1(15) 

DO 10 J«l,15 
KAY-KAG(J) 

F1=CALF1(J) 

TLC=2.*L/C 

TLCK=(1.339*TLC*KAY)**2 

PHISIG=(TLC**(5./3.) ) *{l.+(8./3.)*TLCK)/( (1.+TLCK)**(11./6.)) 
F(J)»(KAY**2) *F1*PHISIG 
10 CONTINUE 
RETURN 
END 

SUBROUTINE CALC 
IMPLICIT REAL*8(A-H,K,0-Z) 

REAL*8 HU,H,NRSQ,NISQ,L 
INTEGER APM 



COMPLEX* 16 DCMPLX , DCONJG , YDOT , YDTC , Y , YTOP , YBOT 

COMPLEX* 16 AlCOM , A2COM , BICOM , B2COM , CICOM , C2COM 

COMPLEX* 16 SOP , ABPHZ , ABY , PN , PD , YN , YD , DM, DD , PAl , PBl , PA2 , PB2 

COMMON/ INP/H , RHO , VS , AO , ETA , FDEL , PDLDT , W, C , CT , ET , S , ST , AW , AT , PC 

I, RAX(60) , KAPPA, KAY, M,L 
COMMON/FLT/APM , NFC , NHT 

COMMON/ CALCU/V, VC, PI , APRI , BPRI , SMLA, SMLAT, R0 , R1 , E,COl ,C02 ,C03 , R2 , S 

II, MU 

COMMON/ INT/BIGK , AS IG, FI , KPRIME , KC , KO , NO 
COMMON/PRT/YDOTSQ , P2 ,F3 , IP, JP 
COMHON/FAST/TVC,GKl ,GK2,GK3 ,GK4 ,GKPHI,GKY 

COMMON/GEN/TPT(4,10) ,TPB(4,10) ,TYT(4,10) ,TYB(4,10) ,COND(4) ,CONY(4) 
l,CONP(4) ,QUADA(4) ,QUAOB(4) ,TDB(4) ,QUAFA(4) ,QOAPB(4) ,QUAYA(4) ,QDAYB 
2(4) ,QUAPTA(4) ,QUAPTB(4) ,QUAYTA(4) ,QOAYTB(4) ,NS(4) ,MS(4) 
COMMOH/GEN2/K3(4) ,KPHI(4) ,HPT,NPB,NYT,NYB 
IF(KAY.EQ.0.)F1«0.0D0 
IF(KAY.EQ.0.)GO TO 800 
OMEGA«TVC*KAY 
SOP<° DCMPLX ( 0 . 0D0 , OMEGA) 

DN-DCMPLX ( 1 . 0D0 , 0 . 0D0 ) 

DD-DCMPLX ( 0 . 0D0 , 0 . 0D0) 

PN- DCMPLX (1.0D0,0.0D0) 

IF ( APM. EQ. 0 ) PN-DCHPLX ( 0 . 0D0 , 0 . 0D0) 

PD-DCMPLX ( 1 . 0D0 , 0 . 0D0) 

YN- DCMPLX ( 1 . 0D0 , 0 . 0D0 ) 

IF (APM. EQ. 0) YN« DCMPLX ( 0 . 0D0 , 0 . 0D0) 

YD-DCMPLX ( 1 . 0D0 , 0 . 0D0 ) 

IF(APM.EQ.0)GO TO 50 
DO 10 J-1,NPT 

10 PN-PN*(TPT(APM,J)*SOP+l.) 

DO 20 J-1,NPB 

20 PD-PD*(TPB(APM,J)*SOP+l.) 

DO 30 J»1,NYT 
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30 YN=YN*(TYT(APM,J)*S0P+1.) 

DO 40 J»1,NYB 

40 YD»YD*(TYB(APM, J) *SOP+l.) 

DD«COND(APH) *GK1*GK2*GK3* (QOADA(APM) *SOP**2+QUADB ( APH) *SOP+l.) 
DN=TDB(APM) ‘SOP/GK4+1. 

PN=CONP(APH) *GKPHI* (QUAPTA(APM) *SOP**2+QUAPTB ( APM) *SOP+l.) *PN 
PD= (QUAPA ( APM) *SOP** 2+QUAPB ( APM) *SOP+l . ) * ( SOP**NS ( APM) ) *PD 
YN=CONY ( APM) *GKY* ( V/57 .3 ) * (QUAYTA(APM) *SOP**2+QOAYTB ( APM) *SOP+l . ) * 
lYN 

YD*(QUAYA(APM) *SOP**2+QUAYB(APM) *SOP+l .) * (SOP**MS (APM) ) *YD 
50 STAB»FDLDT*C02*KAY 

TAIL—DCMPLX (FDEL , STAB) 

ABPHI-TAIL* (DD/DN) * (PN/PD) 

ABY-TAIL* (DD/DN) * (YN/YD) 

PARI «=- 2 . * MU * TVC * ( K AY * * 2 . ) 

PAIl-TVC* (R0+R1* ( 1 ,-ETA) ) *KAY 
PBRl— (R0+R1* (l.-ETA) ) 

PBI1=(R0*(2.*CO3-1.)-R1*R2)*KAY 

C1«R0+R1*(DCOS(KAY*S1)-ETA) 

CIBAR— R1*DSIN(KAY*S1) 

PAR2-0.0D0 

. PAI2-TVC* (CO3*R0-Cbl*Rl* (l.-ETA) ) *KAY 
PBR2— (CO3*R0-CO1*R1* (l.-ETA) +4 . *WJ* (KAPPA**2) * (KAY**2)/(C**2) ) 
PBI2»(CO3*R0*(2.*CO3-1.)+CO1*R1*R2)*KAY 
C 2 *C03 * R0 -COl * R1 * ( DCOS ( K AY* SI ) -ETA) 

C2BAR=C01*R1*DSIN(KAY*S1) 

PAl =DCMPLX ( PARI , PAIl ) 

PB1=DCMPLX(PBR1,PBI1) 

PA2=DCMPLX (PAR2 ,PAI2) 

PB2=DCMPLX(PBR2,PBI2) 

A1C0M=R1*ABY+PA1 

A2C0M«-C01*R1*ABY+PA2 


B1C0M=R1*ABPHI+PB1 



B2C0H—C01*R1*ABPHI+PB2 
CICOM-DCMPLX (Cl ,C1BAR) 

C2COM»DCMPLX (C2 ,C2BAR) 
yTOP»ClCOM*B2COM-C2COM*BlCOM 
YBOT-A1COM*B2COM-A2COM*B1COM 
y-YTOP/YBOT 
700 YDOO>TVC*KAY*Y 
YDTC-DCONJG(YDOT) 

YDOTSQ-YDOT*YDTC 

BETASQ-1.-(M**2) 

F2»1./(1.+(1.5*PI*KAY/BETASQ)+((PI**2)*M*(KAY**2)/BBTASQ) ) 

F3-1 ./ ( 1 .+ ( . 55 *AW*KAY) ) 

Fl-4 . * (MU**2) * (KAY**2) *YDOTSQ*F2*P3 
800 RETURN 
END 

SUBROUTINE PRINT 
IMPLICIT REAL*8(A-H,K,0-Z) 

REAL*8 MU,M,L,NO 
INTEGER APM 

COMMON/ INP/H , RHO , VS , AO , ETA , FDEL , FDLDT , W, C , CT , ET , S , ST , AW, AT , FC 

I, KAX(60) , KAPPA, KAY, M,L 
COMMON/FLT/APM , NFC , NWT 

COMMON/CALCU/ V, VC , FI , APRI , BFRI , SMLA, SHLAT , R0 , R1 , E , COl , C02 , C03 , R2 , S 

II, MU 

COMMON/ INT/BIGK , AS IG , F 1 , KFRIME , KC , KO , NO 
COMMON/PRT/YDOTSQ, F2 ,F3 , IP, JP 
COMMON/FAST/TVC , GKl ,GK2 ,GK3 ,GK4 ,GRPHI ,GKY 
DIMENSION F1X(60) ,YDTSQ(60) ,F2X(60) ,F3X(60) ,FISIG(60) 

IPRT-0 

IF (NWT. NE. 2) GO TO 20 
IPRT-1 

DO 10 J»l,41 
KAY=KAX(J) 
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CALL CALC 


F1X(J)=F1 

TLC=2.*L/C 

TLCK-(1.339*TLC*KAY)**2 
YDTSQ(J) =YDOTSQ 
F2X(J)»F2 
F3X(J)»F3 

FISIG(J) »(TLC**(5./3.) ) *(l.+(8./3.) *TLCK)/( (l.+TLCK)** (11 ./6.)) 
10 CONTIKUE 
20 CONTINUE 
JP-JP+1 
KP=(JP-l)/3 
KP=3*KP+1 

IF(JP.EQ.KP)IP-IP+1 
IF(JP.EQ.KP)WRITE(6,30)IP 
CALL CODE (NFC, NHT,APN) 

AS IG- VC *B IGK/ ( 3 2 . 2 *MU ) 

ASIGVC-ASIG/VC 

WRITE(6,40)M,H,V,VC,MO,KC 

WRITE(6,50)AO,SMLA,SMLAT 

WRITE(6,60)R0,R1,R2,S1 

WRITE(6,70) BIGK,KPRIME,ASIG,ASIGVC,KO,NO 

IF(IPRT.EQ.0)GO TO 100 

IP=IP+1 

JP»0 

WRITE (6, 30) IP 

CALL CODE ( NFC, NWTjAPM) 

WRITE (6, 80) 

DO 90 J«l,41 
KAY>=KAX(J) 

F2F3=F2X(J) *F3X(J) 

FUNF1=F1X(J)*FISIG(J) 

FUN1F1»(KAY**2) *FUNF1 
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GFPSD-KAY*F1X(J)*PISIG(J) 

GNFPSD- ( ( VC/ ( 3 2 . 2 * MU ) ) * * 2 ) *G PPSD/ ( PI * (TLC* * ( 2 . /3 . ) ) ) 

ANFPSD-GFPSD/ ( BIGK**2) 

WRITE(6,85)KAY,YDTSQ{J) ,F2F3,F1X(J) ,FISIG(J) ,FUNF1,FUN1F1 ,GFPSD,GN 
iFPSDjAMFPSD 
90 CONTINUE 
100 RETURN 

30 FORMAT ( '1' ,10X, 'GUST RESPONSE INTEGRATIONS ' ,10X PAGE' , 13///) 

40 FORMAT(6X, 'M' ,10X, 'H' ,10X, 'V ,9X, ' V/C »8X, ' MU ' ,9X, 'KC' ,/,lX,lP6Ell 
1.3/) 

50 FORMAT(5X, 'AO' ,10X, 'A' ,9X, 'AT' ,/ ,1X,1P3E11 .3/) 

60 FORMAT (5X, 'R0' ,9X, 'R1 ' ,9X, 'R2 ' ,9X, 'SI ' ,/,lX,lP4Ell .3/) 

70 FORMAT(6X, 'K' ,8X, 'KPRIME' ,7X, 'ASIG' ,5X, 'ASIG/V/C ,6X, 'KO' ,9X, 'W' » 
1/,1X,1P6E11.3//////) 

80 FORMAT (5X, 'KAY' ,8X, 'YDOTSQ' ,5X, 'F2*F3' ,7X, 'FI' ,7X, 'PHI/SIG' ,2X, 'FI 
1*PHI/SIG' ,2X, 'KSQ*F1*PHI/SIG* ,3X, 'GFPSD* ,5X, 'GNFPSD' ,5X, 'ANFPSD' ,/ 
1) 

85 F0RMAT(1X,1P6E11.3,3X,1PE11.3,2X,1F3E11.3/) 

. END 

SUBROUTINE CODE (NFC, NWT,APM) 

INTEGER APM 
- WRITE (6, 800) 

IF(NWT.NE.l)GOTO 620 . 

WRITE (6, 600) NFC, APM 
GO TO 700 

620 IF(NWT.NE.2)GOTO 660 
WRITE (6, 640) NFC, APM 
GOTO 700 

660 WRITE (6, 6 80) NFC, APM 
700 WRITE(6,800) 

RETURN 

600 FORMAT (lx, 'CASE NUMBER' ,2X,I2 , '-A-' ,11) 

640 FORMAT (IX, 'CASE NUMBER' ,2X, 12 , '-B-' ,11) 




680 FORMAT (IX, 'CASE NUMBER' ,2X, 12 
800 FORMATdX,' 
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0.84 

30 

40 
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28 

32 

0.950 
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0.749 

0.741 

0.741 
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0.719 

0,721 
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0.118 

0.110 
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0,094 
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0.111 

0.090 

0.120 
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0.030 

0.030 

0,157 

0il23 

0.096 
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0.101 

0 

0.274 

0.279 

0.239 

0.106 

0.098 

0.112 

0.082 
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0.108 
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Table 1 Airplane Characteristics 
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7.- Generalized AgAPF versus mass ratio parameter, A/P mode 
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Figure 8.- Generalized A^APF versus mass ratio parameter, A/P mode Cruise. 
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Figure 9.- Generalized AqAPF versus mass ratio parameter, A/P mode Turbulent 
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Figure 12.- Acceleration increment zero crossing number versus mass ratio 

parameter, A/P mode Off, V/c = 40 sec“^. 
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o Figure 16.- Zero crossing nvimber APF versus mass ratio parameter, A/P 

mode Cruise, V/c = 40 sec~^. 
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Figure:. 19..- GNFPSD versus frequency ratio, Cliinb ami cTescent, A/P mode Off 
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Figure 20.- GNFPSD versus frequency ratio. Climb and descent, A/P mode Climb/Descent. 
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Figure- 2T.- GNFPSD versus frequency ratio, Eow cruise, S/P' mode Off . 
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Figure 22.- GNFPSD versus frequency ratio. Low cruise, A/P mode Cruise. 
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'Fiqure 23-. - t3NFPSD versus frequency -ratio , Low cruise , 7L/P mode "Turbulent . 
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o. High cruise, A/P mode Off 
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Figure 25,- GNFPSD versus frequency 
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Figure 26.- GNFPSD versus frequency ratio. High cruise, A/P mode Turbulent. 
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Figure -27. - -ANFPSD versus frequency oratio. Climb and -descent, A/P mode Off 
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Figure 28.- ANFPSD versus frequency ratio. Climb and descent, A/P mode Climb/Descent. 
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Figure 30.- ANFPSD versus frequency ratio. Low cruise, A/P mode Cruise. 
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Fiqurje 31,- •Al'JFJPSD versus frequency ratio. Low crnise, A/P mode Turbulent. 
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Figvire 32.- ANFPSD versus frequency ratio. High cruise, A/P mode Off. 
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Pigure 33.- ANFPSD versus frequency ratio. High cruise, A/P mode Cruise. 
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Figure 34.- ANFPSD versus frequency ratio. High cruise, A/P mode Turbulent . 
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